Feasibility of the finite amplitude method in covariant density functional theory 
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Self-consistent relativistic random-phase approximation (RPA) in the radial coordinate repre- 
sentation is established by using the finite amplitude method (FAM). Taking the isoscalar giant 
monopole resonance in spherical nuclei as example, the feasibility of the FAM for the covariant den- 
sity functionals is demonstrated, and the newly developed methods are verified by the conventional 
RPA calculations. In the present relativistic RPA calculations, the effects of the Dirac sea can be 
automatically taken into account in the coordinate-space representation. The rearrangement terms 
due to the density-dependent couplings can be implicitly calculated without extra computational 
costs in both iterative and matrix FAM schemes. 
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I. INTRODUCTION 

By reducing the quantum mechanical many-body 
problems formulated in terms of A^-body wave functions 
to the one-body local density distributions, the density 
functional theory (DFT) of Kohn and Sham [l| has ac- 
complished great success in many different fields of mod- 
ern physics. No other method achieves comparable accu- 
racy at the same computational costs. In nuclear physics, 
the DFT has been widely used since the 1970s [2|. In par- 
ticular, its covariant version in the relativistic framework 
has received much attention during the past decades. 

The covariant density functional theory (CDFT) 0, Q 
takes the Lorentz invariance into account. In this frame- 
work, the representation with large scalar and vector 
fields, of a few hundred MeV, provides a consistent treat- 
ment of the spin degrees of freedom. The Lorentz sym- 
metry leads to the unification of the time-even and time- 
odd components in the corresponding functionals. The 
Lorentz symmetry also puts stringent restrictions on the 
number of parameters without reducing the quality of 
the agreement with experimental data. Over the years, a 
large variety of nuclear phenomena have been described 
successfully by the CDFT [5|-|3|, including the equation 
of state in symmetric and asymmetric nuclear matter, 
ground-state properties of finite spherical and deformed 
nuclei all over the nuclear chart, collective rotational 
and vibrational excitations, fission landscapes, low-lying 
spectra of transitional nuclei involving quantum phase 
transitions in finite nuclear systems, and so on. 

Focusing on the vibrational excitations, the random- 
phase approximation (RPA) Q is one of the leading the- 
ories applicable to both low-lying excited states and gi- 
ant resonances. In the relativistic framework, the self- 
consistent and quantitative RPA calculations were real- 
ized after recognizing the importance of the Dirac sea 



[l^-Uyl- It has been proved that the relativistic RPA 
is equivalent to the corresponding time-dependent rela- 
tivistic mean-field (RMF) theory in the small amplitude 
limit, only if the particle-hole (ph) configurations include 
not only the pairs formed from the occupied and unoc- 
cupied Fermi states but also the pairs formed from the 
Dirac states and occupied Fermi states |12j . 

From then on, great efforts have been dedicated to de- 
veloping the self-consistent RPA approaches in the rel- 
ativistic framework [7|. The formalism for the nonlin- 
ear meson-exchange interactions can be traced back to 
Refs. [13, [l3 ■ For the density-dependent meson- nucleon 
couplings, the explicit rearrangement terms in the ph 
residual interactions have been derived |19j . The rel- 
ativistic quasiparticle RPA (QRPA) [201 ^^ been de- 
veloped based on the canonical single-nucleon basis of 
the relativistic Hartree-Bogoliubov theory for giant res- 
onances [2ll423l |. pygmy resonances j24| - |26| . and low- 
lying vibrational states ||3-[2^. The relativistic (Q)RPA 
has also been extended to the charge-exchange chan- 
nels [30, |311 for the nuclear spin-isospin resonances p32| - 
[3a | , /3-decay rates [33, [S^ , muon-capture rates [SO] , and 
neutrino-nucleus reactions [io, |41| . In addition, the rel- 
ativistic RPA with finite temperature [42|, \^ and the 
continuum (Q)RPA [Jj - liTJ have been estabhshed. To 
go beyond the mean field, the particle vibrational cou- 
pling has also been taken into account [i^, H^ ■ 

Recently, a fully self-consistent relativistic RPA [50| 
has been established based on the relativistic Hartree- 
Fock theory [5ll - l53| . It is shown that not only the 
Gamow- Teller resonances but also the fine structure 
of spin-dipole resonances can be well reproduced with- 
out any readjustment of the energy functional [50J, |5J]. 
This self-consistent RPA has also been applied to eval- 
uate the isospin symmetry-breaking corrections to the 
supcrallowed /3 transitions for the unitarity test of 



Cabibbo-Kobayashi-Maskawa matrix [55|- The corre- 
sponding QRPA [5Q] based on the relativistic Hartree- 
Fock-Bogohubov theory [S^ has been developed and 
used for a systematic study of /3-decay half-hves of 
neutron-rich even-even nuclei with 20 ^ Z ^ 50, where 
the isospin-depcndent isoscalar proton-neutron pairing is 
found to play a very important role. 

However, the above investigations are essentially re- 
stricted within the spherical symmetry. The conventional 
RPA calculations in the matrix form face a big compu- 
tational challenge when the number of ph configurations 
Nph becomes huge as in the deformed cases. So far, the 
only self-consistent deformed (Q)RPA in the relativistic 
framework was developed by Pena Arteaga et al. |58l.l59|. 
Note that, even in the non- relativistic framework, the de- 
formed (Q)RPA in the matrix form is also a hard task. 
There are a few recent attempts for the Skyrme ene rgy 
density functionals in the axially symmetric case |60l464l | 
and in the triaxial case [65|, as well as for the Gogny en- 
ergy density functionals in the axial case [6a]. The full 
three-dimensional calculations have been carried out only 
using the real-time methods J67H7C| . 

As a promising solution for this computational chal- 
lenge, the so-called finite amplitude method (FAM) was 
proposed in Rcf. [7l[. In this method, the effects of 
residual interactions are evaluated in a numerical way by 
considering a finite density deviation around the ground 
state. In such a way, the self-consistent RPA calcula- 
tions become possible with a little extension of the static 
Hartree(-Fock) code. Furthermore, by using the iter- 
ative methods for the RPA equation, the computation 
time is close to a linear dependence on Nph, instead of 
a dependence between N^f^ and A^^^ in the diagonaliza- 
tion scheme [T^l- This advantage is crucial when Nph 
becomes huge. In the non-rclativistic framework with 
Skyrme energy density functionals, the feasibility, accu- 
racy, and efficiency of FAM have been demonstrated for 
the RPA in the three-dimensionally deformed cases in the 
coordinate-space representation [7l|, [7^, UM ^^^ for the 
QRPA in the spherical (zl, ^ and axially deformed [zi] 
cases in the quasiparticle-basis representation. Iterative 
algorithms for (Q)RPA solutions have also been devel- 
oped recently, based on the Arnold! process |77l . [78j and 
on the conjugate gradient method |79j . The readers arc 
also referred to Ref. [80] for a recent review. 

Therefore, it is worthwhile to develop the self- 
consistent relativistic RPA by using the finite amplitude 
method. In particular, special attentions should be paid 
to the unique features of covariant density functionals, in- 
cluding the effects of the Dirac sea and the rearrangement 
terms for the density-dependent interactions. These re- 
arrangement terms are usually more sophisticated than 
those in the Skyrme functionals, and cause heavy com- 
putations [13 • On the other hand, the covariant density 
functionals hold the Lorentz invariance, which leads to 
the unification of their time-even and time-odd compo- 
nents. This makes the modification in the ground-state 
code straightforward. 



In this work, our premier purpose is to verify the fea- 
sibility of the FAM in the CDFT, with special attentions 
to the Dirac sea and the rearrangement terms. For a ba- 
sic demonstration, the self-consistent RPA is established 
based on the spherical density-dependent point-coupling 
RMF theory by using the FAM. 

The paper is organized as follows: In Sec. |lll the key 
formulas of the density-dependent point-coupling RMF 
theory and the corresponding self-consistent RPA, and 
the formalism of both iterative and matrix FAM will be 
presented. In Sec. lIIIl the numerical details will be shown 
with the main focus on the boundary conditions of the X 
and Y amplitudes in the coordinate-space representation. 
In Scc. lIVi a benchmark test will be given and the effects 
of the box size, Dirac sea, and rearrangement terms on 
the isoscalar giant monopole resonances (ISGMR) will be 
discussed. Finally, a summary will be given in Sec. |Vl 



II. THEORETICAL FRAMEWORK 
A. Point-coupling relativistic mean-field theory 

Successful CDFT can be traced back to the RMF mod- 
els introduced by Walecka and Serot Q- Since then, the 
popular RMF models [3-|a| are based on the finite-range 
meson-exchange representation, in which the nucleus is 
described as a system of Dirac nucleons that interact with 
each other via the exchange of mesons. 

Recently, the CDFT framework has been reinterpreted 
by the relativistic Kohn-Sham scheme, and the function- 
als have been developed based on the zero-range point- 
coupling interactions [8l|. In this framework, the meson 
exchange in each channel is replaced by the correspond- 
ing local four-point contact interaction between nucleons. 
Such point-coupling model has attracted more and more 
attentions during the pastyears due to its simplicity and 
several other advantages J8l l8^490l |. In particular, for the 
present study, by directly expressing the mean-field po- 
tentials in terms of nuclcon densities and currents, the 
FAM can be applied in a more straightforward way. 

In this section, we recapitulate the key formulas of the 
point-coupling RMF theory for the FAM calculations, 
in particular, those related to the currents and space- 
component of the Coulomb field. 

The effective Lagrangian density of the density- 
dependent point-coupling RMF theory reads [8^ 



ei^^>^A,^^^i^-\F^'^F,^, 



(1) 



where M is the nucleon mass, and the field tensor for 
photons reads Fi^" = di^A" - d^Af". While the couphng 
parameter 6s is a constant, the coupling strengths of the 



four-nucleon interactions in the scalar (S), vector (V), 
and isovector-vector (tV) channels are analytical func- 
tions with respect to the baryonic density pf,, 



as{pb) ^ as + ibs + csx)e '^^'^ , 
aviPb) ^ av + bye''^^'', 
atviPb) = 6tve-'^'''^ 



(2a) 
(2b) 
(2c) 



with X — Pb/ Psnt, and psat denotes the saturation density 
of symmetric nuclear matter. 

In this paper, the vectors in coordinate space are de- 
noted by bold type, and vectors in isospin space are 
denoted by arrows. Greek indices p, v run over the 
Minkowski indices 0, 1, 2, and 3. 

The effective Hamiltonian H can be obtained with the 
general Legendrc transformation. Together with the trial 
ground state |<f>o) a-s a Slater determinant, as well as 
the Hartree and no-sea approximations, the energy func- 
tional can be written as 

E = ($o| H |$o) = E,, + Es + Ev + Etv + Ea, (3) 

where the first term is the kinetic energy, and the oth- 
ers correspond to contributions from the scalar, vec- 
tor, isovector-vector channels and Coulomb field, respec- 
tively. 

The Dirac equation for nucleons, 



together with the additional rearrangement term due to 
the density-dependent coupling strengths, 



^i? = t: 



das 2 
opb 



da 



dpi 



V .,1 . 



da 



dp, 



tv -^n -f 
Jtv ■ JtVp, 



(6) 

In the above expressions, ps, jy, j^y, and A'^ are the 
scalar density, the isoscalar and isovector four-currents, 
as well as the Coulomb field, respectively. The nuclear 
baryonic density pb corresponds to the time-component 
of the isoscalar four-current jy . 

It is worthwhile to emphasize here that the space- 
components of the four-currents and the Coulomb field 
must be kept explicitly for the following applications of 
FAM, even though they in general vanish in the ground 
state of systems with the time-reversal symmetry, e.g., 
even-even nuclei. 

For the systems with spherical symmetry, the single- 
particle wave functions have the form of 



V'o 



^Gair) 
Fa{r)a ■ 



^air)xi{qa) 



(7) 



h\lpa) 



IV'a) 



(4) 



is then obtained by the variation principle. The one-body 
mean-field Hamiltonian h is composed of the kinetic term 
/ifc, the scalar ^15, vector hy, isovector-vector hty, and 
Coulomb Ha terms, i.e.. 



h^ = ~ia ■ V + 7°Af , 
hs = j°{asps + <5sA/9s), 

hy = 'y°j^ayj!^, 



htv 
Ha 



l°lt.atvT-j^y, 



e7°7M 



;i-r3) 



A". 



(5a) 
(5b) 
(5c) 
(5d) 

(5e) 



where ^ "„ (f ) arc the spherical harmonics spinors, 
Xi('Za) the isospinors. The single-particle eigenstates are 
specified by the set of quantum numbers a = {a, ma) = 
iqa,'n'a,la,ja,'m'a), and the good quantum number Ka = 
T(ja + 1/2) for ja = ^a ± 1/2. Within this phase conven- 
tion between the upper and lower components, the wave 
functions G{r) and F{r) can be simultaneously chosen 
as real functions for the ground-state descriptions. In 
contrast, for the FAM built beyond, both G{r) and F{r) 
become complex functions, so one should be careful to 
distinguish them from their complex conjugates G*{r) 
and F*{r) from the very beginning. 

The radial Dirac equation reads 






(8) 



with the scalar and vector potentials 



^sir) = asPs{r)+Ss (ps(r) + -p's{r)] , 

1 - T3 
Eo(r) = aypy{r) + atyptv{r)T3 + e — - — Ao{r) + Sfl(r), 

1 - T3 
Sy(r) = ayjy{r) + atvjtv{r)T3 + e — - — Ay{r). 



(9a) 
(9b) 
(9c) 



The rearrangement terms only contribute to the time-component of the vector potential, which read 



^-<-'4{|f^^« + ^<^v.M+-'-^'-)) 



da 



tV , 2 



dpb 



ipiv{r)+jtvir)) 



(10) 



r 



The densities and currents are expressed as 
p'-t^ -J^ll^l \Gl{r)Ga{r) - KM^Jr)] , (11a) 



p'r 



J^ 



\ E3a {G:Xr)Ga{r) + F:(r)F^(r)\ , (lib) 



47]- 7-2 



475-7-2 



Y.fAGl{r)F,{r)-F:{r)Ga(r)\^ (lie) 



as a linear response to the weak perturbation. In the 
frequency representation, the above equation is expressed 
as 



uj^piuj) = [/io,<5/9(w)] + {^KlS) + Kxt(w),/Oo]- (15) 
In practical calculations, it is convenient to adopt 
the single-particle (Kohn-Sham) orbitals to represent the 
density matrix. 



with 2\ = 2ja + 1- The isoscalar densities and currents 
are the sum of the neutron and proton contributions, 
while the isovector ones are the differences between the 
neutron and proton contributions. The Coulomb fields 
arc calculated with the Green's function method, i.e.. 



Ao(r) = e/drV'pl?^(r')- 



'r> 
where r> = max{r, r'} and r< = min{r, r'}. 



(12a) 
(12b) 



p(t)=5^|^,(t))(V^,(t)| 



(16) 



As a result, the density deviation in the frequency repre- 
sentation can be expressed as 



6p{co) ^ E{|X,(..)) (0,1 + 10,) {Y,{uj)\} (17) 



B. Linear response and random-phase 
approximation 

The RPA equation is known to be equivalent to the 
time-dependent Hartree(-Fock) equation in the small am- 
plitude limit [9[. In order to make the FAM clear in 
the next section, we first briefly recall the derivation of 
the standard RPA equation by following the notations in 
Ref. [nl- 

The static Hartree or Hartree-Fock equation. 



with the so-called forward X{u!) and backward Y{u!) am- 
plitudes and the occupied eigenstates {|0i)} of hg in 
Eq. ^. It is slightly tricky that one must take the ket 
\Xi{uj)) and bra (yi(w)| states independent, since Sp{uj) 
is not Hermitian. But, this point is in fact well known as 
the solutions of the RPA equation shown below. Here- 
after, \(j)a) represent the eigenstates of ho, and indices i,j 
(to, n) run over the hole (particle) states. 

By expanding the X(u;) and Y{uj) amplitudes on the 
basis of particle states. 



[h[p],p]=0, 



(13) 



determines the ground-state density p = po satisfying 
p^ = p, and the one-body mean-field Hamiltonian ho ~ 
h[po]. 

When a time-dependent external perturbation Vc-xt{t) 
is present, the density deviation Sp{t) = p{t) — pQ obeys 



i-Sp{t) = [ho,6p{t)] + [5h{t) + V,^t{t),Po] 
at 



(14) 



m>A 



(18a) 
(18b) 



m>A 



one can derive the well-known RPA equation in the ma- 
trix form, 



B* A* 

mi.nj "^mi^nj 



1 
-1 



71 J \^) 



X. 






(19) 



The RPA matrices A and B and vectors / and g read 

dh 



A. 



S. 



■rni^nj 

'7ni,nj — 
J mi 



(e„i - ei)S„inSij + (0m I 



dp. 



'nj 



I't'i) = (e™ - <^i)SmnSij + {4'm4'j\ Vph |0„(/'i) , 



P=Po 



dh 



\4'i) = (^rn'/'nl Vph \4>]4'i) , 



dpjn 

Kxt(w) \4>i) , gmi = (0j| 14xt(w) |0m) . 



P=Po 



(20a) 

(20b) 
(20c) 



For the self-consistent RPA calculations [9|, the particle-hole residual interactions Vph should be strictly derived 
from the second derivative of the energy functional E shown in Eq. ([3l. The ph residual interactions for the point- 
coupling RMF theory with nonlinear couplings can be found in Ref. [9l[. In contrast, the density dependence in the 
coupling strengths a introduces additional rearrangement terms in Vph uM- Explicitly, the ph residual interactions 
are composed of 



F/,(l,2) = |a5[7"]i[7°]2 + ^P5([7°]i[I]2 + [I]i[7"]2) + 
F^^(l,2) = |av[7V]i[7%M]2 + 2^Pv[I]i[I]2 
F,V^(1,2) = (a*y[7Vr]i • [I'l^^h + ^PtviirzWh 



2^pf 



plmi]2-Ssh°S7]i ■ [7°V]2 U(ri -r2), (21a 



1 d^av 
2^ 



p2,[I]i[I]2U(ri-r2), 



1 



^A(^.'^-rJ-'°-''-'-r^'^'°''--'-r^'W^i 



Jl[^3j2j 



1 d'^atv 

2 dpi 



p2^[I]i[I]2U(ri-r2), 



(21b) 
(21c) 
(21d) 



where I denotes the 4x4 unit matrix. The rearrange- 
ment terms correspond to those containing da/dpb or 
d^a/dpl- They are calculated term by term separately 
in the conventional RPA calculations. 

Meanwhile, it is also important to emphasize the effects 
of the Dirac sea. The relativistic RPA is equivalent to 
the time-dependent RMF theory in the small amplitude 
limit, only when the particle states to, n include not only 
the states above the Fermi surface but also the states in 
the Dirac sea [l^l ■ It is due to the no-sea approximation 
used in the ground-state calculations. In other words, the 
ensemble of all these unoccupied states together provides 
a complete set of basis for particle states. 



C. Iterative finite amplitude method 

In Ref. [7l|, the FAM was proposed as a simpler and 
more efEcient approach to the solutions of the linear re- 
sponse equation (|15p . This method docs not require ex- 
plicit evaluation of the residual interactions 5h/ 5p as in 
Eq. (PO)) . Instead, by multiplying with the ket \(j)i) and 
bra {4>i\ of only hole states on both sides of Eq. ([T5|). 
respectively, one has 

w \X,{io)) = {ho - eO \X,{uj)) + Q(14xt(^) + 5h{Lu)) |0,;) , 

(22a) 

u* \Y,{u)) ^ -{ho - e,;) \Y,{u)) - Q{V,U^) + Sh\uj)) 10,) 

(22b) 



where Q = 1 — X); I0j) (0il i^ ^ projection operator onto 
the particle space. 

The induced fields Sh{uj) and Sh^{uj) shown above are 
calculated by using the following finite difference with a 
sufficiently small number rj: 



sh{u;)^-{h[{^i;'\,m-hm,m) (23) 



with (?A^| = 
and 



+ Tj{Yi{io)\ and \^,) = 10,) +77|X,(a;)), 

Sh^c.) = -{h[{y^'\,m~h[{<j>\,m) (24) 

V 

with (V','! = (0.1 + V {M^)\ and IV'.) = 10.) + V \Y^{uJ)). 

For the present calculations with spherical symmetry, 
it is convenient to rewrite Eq. (|22p in coordinate space. 
Assuming the monopole perturbation, 

Voxt{r,uj) = Vcy,t{r,uj)Yoo{r), (25a) 

X,{r) ^ X,{r)Yoo{r), (25b) 

y.(r) = y.(r)yoo(f), (25c) 

the corresponding radial FAM equations read 

Q [{ho{r) - e, - oj)Xi{r,uj) + 6h{r, cj)0i(r)] 
= -QKxt(r,w)0.(r), (26a) 

Q [{ho{r) - e, + UJ*)Y^{r,iu) + Sh\r,iu)cf>,{r)]* 



vUr,u)Ur) 



(26b) 



In the relativistic framework, h{r) is a 2 x 2 matrix as 
shown in the radial Dirac equation (|S]), and (/>i(r) = 
{Gi{r) Fi{r))'^ . Therefore, the X and Y amphtudes are 
also composed of the upper and lower components. 



Mr) 



Xcrir) 
XF^{r) 



Y^ir) 



Ycrir) 
YF,{r) 



(27) 



As emphasized in the previous section, the effects of 
the Dirac sea must be taken into account, which is ex- 
pressed in an explicit way in the conventional expansions 
(fT8| . In contrast, here the X and Y amplitudes are ex- 
panded on the mesh points {rk} in coordinate space. In 
such a way, on one hand, the effects of the Dirac sea can- 
not be identified or isolated; on the other hand, from 
the mathematical point of view, the coordinate space 
J2r |i") ("^1 ~ J2j I'f'j) {4'j\^ can also provide a complete set 
of basis for particle states. 

The induced fields 5h{r) and (5/i^(r) are evaluated by 
using Eqs. ([25)1 and (P^ . The procedure in practice is as 
follows: with a given set of {Xi{r)} and {Yi{r)}, one cal- 
culates the nucleon densities and currents [Eq. ([TT]) ]. new 
couphng strengths [Eq. ^], Coulomb fields [Eq. (fT2)) ]. 
rearrangement self-energy [Eq. ([T0|) ]. scalar and vector 
potentials [Eq. ([9|)], and then the one-body Hamiltonian 
h{r) [Eq. ([S])], sequentially. Since now the X{r) and Y{r) 
amplitudes are independent due to the non-Hermitian 
nature of 6h{r) and 6h'^{r), it is clear that the nucleon 
currents are no longer vanishing. This is the reason why 
these time-odd terms must be kept from the beginning. 

In order to include both the normal and rearrangement 
terms in the ph residual interactions as explicitly shown 
in Eq. (PT|) . one simply needs to re-calculate the coupling 
strengths a appearing in Eq. © and their derivatives 
da/dpb in Eq. (fTO)) by using Eq. ([2]) for each given set of 
{Xi{r)} and {y,;(r)}. If one skips this step, i.e., keeps a 
and da/dpb always unchanged, the consequence is that 
the normal terms in Vph remain, but all of the rearrange- 
ment terms are neglected. 

This FAM equation is a standard linear algebraic equa- 
tion of the form. Ax = b, which can be solved within the 
iterative scheme. In such a way, we do not need to con- 
struct the matrix elements of A explicitly, but only to 
evaluate Ax for a given vector x. In the following, we 
denote this iterative finite amplitude method as i-FAM. 

Adopting the w-independent local external field 
Vcxtir,uj) = 0{r), the corresponding transition strengths 
can be calculated with the solutions of Eq. (P5| as 



du> 



-^l^Y^^jf J dr{4{r)0\r)X,{r,uj) 

i 

+ YHr,u;)0\r)Ur)}. (28) 



D. Matrix finite amplitude method 

We introduce another usage of FAM, the so-called ma- 
trix finite amplitude method (m-FAM) shown in Ref. [T^j . 
In this method, the RPA matrices A and B are explicitly 
constructed, but the tedious calculations concerning the 
ph residual interactions Vph in Eqs. ([20]) and (|2ip can be 
avoided. 

First of all, both the occupied and unoccupied eigen- 
states oi ho, {|0i)} and {\4>m)}; are calculated at the 
ground state. Then, instead of dealing with Vph, the ker- 
nels dh/dp in Eq. ([20]) are directly calculated with finite 
difference provided the real parameter 77 is small enough 
to neglect the higher-order terms, i.e., 



dh 



dpn 



= l{h[{i''\,m-hM,\ 

7] 



(29) 



The key point here is to keep all (V'j'l = {(j>i\ and IV'i) = 
\(t>i) unchanged, except for the specific orbital j which 
slightly mixes with another specific orbital n as \iJjj) ~ 
|(/)j) -I- 77 10„). In the same way. 



dh 



dpjn 



P=PO 



1 

V 



{h[{ij'\,m-hm,m 



(30) 



by keeping all {^p^\ = {(j)i\ and \ipi) = \(j)i) unchanged, 
but slightly mixing specific orbitals j with n as (V, I = 
(0, 1 +?7 ((/)„!. 

To include the effects of the Dirac sea, states n run over 
the unoccupied states in both Fermi and Dirac sea. To 
include the effects of the rearrangement terms, one fol- 
lows the same procedure as that in i-FAM shown above. 



III. NUMERICAL DETAILS 

For all the calculations in this paper, the density- 
dependent point-coupling RMF parametrization DD- 
PCl [83| is used and the spherical symmetry is assumed. 
For the ground-state calculations, the radial Dirac equa- 
tion ([H is solved in coordinate space by the fourth- 
order Runge-Kutta method, also known as the shooting 
method, within a spherical box with a box radius R and 
a mesh size dr [92| . The mesh size is fixed as dr = O.I fm, 
while the choice of box size R will be discussed below. 

For the conventional RPA and m-FAM calculations, 
the single-particle energy truncation for constructing the 
RPA matrices A and B in Eq. ^ is [-M, M+200 MeV] , 
i.e., all the bound states in the Dirac sea are taken into 
account. As an example, the corresponding number oi ph 
configurations Nph for J'^ = 0+ excitations in ^osp^-, jg 
1355 with i? = 25 fm, where 524 of them are formed with 
the particle states in the Dirac sea. The convergency 
of this truncation has been examined. The subroutine 
rg . f in EISPACK library is used to diagonalize the non- 
symmetric real RPA matrix. In m-FAM, the parameter 
77 is taken as ry = I0~^. 



For the i-FAM calculations, the frequency uj = E + 
iT/2 contains an imaginary part, and the corresponding 
Lorentzian smearing parameter is F = 1 MeV. The first 
derivative of {Xi{r)} and {Yi{r)} with respective to r is 
performed by the nine-point formula with the boundary 
conditions discussed below. The parameter rj differs for 
every iteration to ensure the linearity |7ll . l73j : 
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m8iii{N{X),N(Y)}' 



NW = J 



\ 



E(^^i^')- 



.(31) 

The truncated version of generalized conjugate residual 
(GCR) method ^, also called ORTHOMIN method, 
is used as the iterative solver, where at maximum 1000 
iterations are stored. The convergent criterion is ||^x — 
^IP/II^IP < 10^^, which provides the relative accuracy 
10"'^ for the transition strengths. 



A. Boundary conditions 

Before further discussions, it is worthwhile to exam- 
ine the boundary conditions of the X and Y amplitudes 
(|27p in the coordinate-space representation. It turns out 
somehow tricky since these amplitudes contain two com- 
ponents instead of one as in the non-relativistic case. 

The boundary conditions for the ground-state radial 
Dirac equation ^ used in the shooting method are 
following [92l: (i) At the origin point, G{r)\r=o = 
F{r)\r=o = 0. (ii) At small distance r — )► 0, G{r) oc 
r~"-,F{r) oc r^~'^ for k < and G{r) ex r^^'^,F{r) ex r** 
for K > 0. (iii) At the box boundary, G{r)\r=R = 0, but 
F{r)\r=R must have a non- vanishing value, otherwise the 
whole wave function will be identically zero. The value of 
F(r)\r=R is eventually determined by the normalization 
condition. 

Accordingly, the consistent boundary conditions of 
X{r) and Y{r) used in i-FAM include: (i) At the ori- 
gin point, XGir)\r=o = ^_F(r)|r=o = >G('')|r=0 = 
Yp{r)\r=o = 0. (ii) At small distance r ^ 0, 
XG{r),YG{r) are odd functions and Xp{r),Yp{r) are 
even functions for even /, while XG{r),YG{r) are even 
functions and Xp{r),Yp{r) are odd functions for odd /. 
(iii) The remaining but critical point is the boundary 
conditions at the box boundary r = R. In addition, out- 
side the box, XG{r)\r>R = Xp{r)\r>R = YG{r)\r>R = 
Yp{r)\ryR, ~ 0, since it is an area out of consideration. 

In order to verify the boundary conditions at r = i?, 
in Fig. [U we show with the solid line the J^ = 0+ unper- 
turbed excitation strengths in ^^O calculated by m-FAM 
with a box radius i? = 20 fni and a mesh size dr = 0.1 fm. 
In m-FAM, the particle states {(j^mir)} correspond to the 
eigenstates of ho{r) with the boundary conditions used 
in the shooting method. Naturally, these boundary con- 
ditions are consistent with the ground-state description. 

In the same figure, the corresponding results cal- 
culated by i-FAM with different boundary conditions 
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FIG. 1: (Color online) The J'^ = 0+ unperturbed excita- 
tion strengths in ^^O calculated by the matrix finite am- 
plitude method (m-FAM) (solid line) with 7? = 20 fm and 
dr — 0.1 fm. The corresponding results calculated by the it- 
erative finite amplitude method (i-FAM) with different box 
boundary conditions are shown with the short-dotted, dash- 
dotted, and dashed lines, respectively. A Lorentzian smearing 
parameter F = 1 MeV is used. 



of X{r) and Y{r) at r = R are shown for compari- 
son. The results obtained by constraining XG{r)\r=R. = 
Xp{r)\r=R = YG{r)\r=R = Yp{r)\r^R = are shown 
with the dash-dotted line, those obtained by constrain- 
ing only XG{r)\r=R = YG{r)\r=R = arc shown with the 
short-dotted line, and those obtained without any con- 
straint at r = R arc shown with the dashed line. The 
tiny difference between the dash-dotted and dashed lines 
shows the effect of changing the box size by one mesh 
point dr, but the visible difference between the short- 
dotted line and the other two is due to the different pre- 
scriptions for the upper and lower components at the 
same position. It can be clearly seen that only the short- 
dotted one with proper boundary conditions is identical 
to the m-FAM result. 

Therefore, the consistent boundary conditions around 
r = R for the X{r) and Y(r) amplitudes in the i-FAM 
calculations read 



XGir)\r^R = Xpir)\r>R = YGir)\r'^R = Ypir)\r 



>R 



(32) 



IV. RESULTS AND DISCUSSION 



In the following discussions, we take the stable and 
radioactive neutron-rich doubly magic nuclei, ^'^^Pb and 
^■^^Sn, as examples. It has been shown that the RMF 
theory can in general nicely reproduce the corresponding 
ground-state properties (see e.g., Ref. |94|). 
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FIG. 2: (Color online) Isoscalar giant monopole resonance 
(ISGMR) in ^°*Pb calculated by m-FAM (solid line) and con- 
ventional RPA (short-dash-dotted line). A Lorentzian smear- 
ing parameter F = 1 MeV is used. 



A. Benchmark tests 

In order to verify the newly developed FAM codes, 
benchmark tests have been performed together with the 
conventional RPA code. The transition strengths of IS- 
GMR in ^°^Pb calculated by m-FAM are compared with 
the conventional RPA results in Fig. [2l where all of the 
common numerical parameters are the same, including 
i? = 25 fm, dr = 0.1 fm, single-particle energy trun- 
cation [-M,M + 200 MeV], and F = 1 MeV. One can 
barely distinguish these two lines in the figure. 

Although we do not show one by one. we have also per- 
formed the conventional RPA calculations for the cases 
without Dirac sea or without rearrangement terms dis- 
cussed below. It is found that all of these results are 
identical to those by the i-FAM and m-FAM calcula- 
tions. This demonstrates the feasibility and accuracy of 
the present FAM codes. 



B. Box size 

In Fig. O the transition strengths of ISGMR in ^ospb 
and ^'^^Sn calculated by m-FAM with box sizes R = 
20, 25, 30, 35 fm are shown with the dashed, solid, dotted, 
and dash-dotted lines, respectively. It is shown that the 
detailed shapes of the resonances change with R to some 
extents. Nevertheless, one of the most important proper- 
ties, the centroid energy mi /mo, docs not depend on R 
up to the digit of 0.01 MeV. Integrating the excitation en- 
ergy from 5 to 25 MeV, the centroid energies in ^°®Pb and 
"2Sn are 14.33 and 16.28 MeV, respectively. The exper- 
imental data in ^o^Pb, mi/mo = 13.96 ± 0.20 MeV [9l, 
can be well reproduced. In the following calculations, the 
box size i? == 25 fm and the mesh size dr = 0.1 fm are 
used. 
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FIG. 3: (Color online) ISGMR in (a) ^'''^Pb and (b) "^Sn cal- 
culated by m-FAM with different box size R. The results cal- 
culated with R — 20, 25, 30, 35 fm are shown with the dashed, 
solid, dotted, and dash-dotted lines, respectively. The experi- 
mental centroid energy in '^"^Pb [93 is denoted by the arrow. 



C. Effects of the Dirac sea 

Comparing with the non-relativistic FAM, it is inter- 
esting to investigate the effects of the Dirac sea in the 
relativistic FAM calculations, in particular, for those us- 
ing the coordinate-space representation. 

The effects of the Dirac sea can be explicitly identi- 
fied in the m-FAM calculations. In Fig. SJ the transition 
strengths of ISGMR in ^''^Pb and ^^^Sn calculated with 
and without the Dirac sea are compared. The results 
including the configurations formed from the occupied 
states in the Fermi sea and unoccupied negative-energy 
states in the Dirac sea are shown with the solid line, while 
the results excluding these configurations are shown with 
the dashed line. It is found that the Dirac sea effects on 
the centroid energies mi /mo of ISGMR in ^^^Pb and 
^^^Sn are as much as 4.00 and 4.26 MeV, respectively. 
This substantial influence is due to the strong couplings 
between the Fermi sea and Dirac sea in the scalar chan- 
nel [13. The experimental data [9a| is reproduced only 
when the Dirac sea is taken into account. 
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FIG. 4: (Color online) ISGMR in (a) ^''^Pb and (b) "^Sn cal- 
culated by i-FAM and m-FAM. The i-FAM results are shown 
with the dotted symbols, while the m-FAM results calculated 
with and without the Dirac sea are shown with the solid and 
dashed lines, respectively. The experimental centroid energy 
in ^"^Pb [95[ is denoted by the arrow. 



FIG. 5: (Color online) ISGMR in (a) ^''^Pb and (b) "^Sn 
calculated by i-FAM and m-FAM. The i-FAM resuhs without 
the rearrangement terms are shown with the dotted symbols, 
while the m-FAM results calculated with and without the re- 
arrangement terms are shown with the solid and dash-dotted 
lines, respectively. 



D. Effects of the rearrangement terms 



In the coordinate-space representation as in i-FAM, 
one can identify no other single-particle eigenstates but 
only the occupied states in the Fermi sea. Just from the 
mathematical point of view, the coordinate space should 
generate another complete set of basis for particle states. 
In Fig. 131 we also plot the corresponding i-FAM results 
with the dotted symbols by taking the energy spacing 
AE = 0.1 MeV. It can be clearly seen that the i-FAM re- 
sults are exactly on top of the m-FAM results that include 
the Dirac sea. This confirms that these two different sets 
of basis are both complete and these two methods are 
equivalent. This also demonstrates that the existence of 
Dirac sea does not introduce additional difficulties for the 
present iterative method in the relativistic scheme, while 
the only price to pay is that the total dimension of the i- 
FAM equations ([26| is now as twice as the non-relativistic 
counterpart. 



It is tedious to calculate the contributions of the re- 
arrangement terms in Vph to the RPA matrix elements 
in the conventional calculations. From Eq. (j2ip . one can 
see that, for one normal term in each channel, there are 
up to 3 rearrangement terms accompanied. In fact, in 
the meson-exchange picture, this number increases to 6 
as shown in Ref. [l^- Even worse, in the RPA based on 
the density-dependent relativistic Hartrce-Fock theory, 
the number of rearrangement terms accompanied can be 
r^ 10^ as a result of an additional summation over the oc- 
cupied orbitals due to the non-locality of the self-energies 

In contrast, as illustrated in Sections III CI and IIIDI 
the effects of the rearrangement terms can be simply 
taken into account in FAM by re-calculating the cou- 
pling strengths a in Eq. (|9]) and their derivatives da/dpb 
in Eq. (|10p with Eq. ([2|) for each given set of {ip'\ and 
\ip). The numerical cost of such a step is totally negligi- 
ble, thus this method is extremely efficient. 
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In Fig. m the transition strengths of ISGMR in ^ospb 
and "'^'^^Sn calculated by m-FAM with and without the 
rearrangement terms arc shown with the solid and dash- 
dotted lines, respectively. Around the main-peak region, 
by taking AE = 0.1 MeV, the i-FAM resuhs calculated 
without the rearrangement terms are also shown with 
the dotted symbols for comparison. The equivalency of 
these two finite amplitude methods is illustrated once 
more, since the rearrangement terms can be switched on 
or off in the same way. Quantitatively, it is found that 
the rearrangement effects on the centroid energies mi/mQ 
of ISGMR in 208pb and ^^asn are 0.53 and 0.26 MeV, 
respectively, which are also substantial. 



V. SUMMARY 

Based on the spherical density-dependent point- 
coupling RMF theory, the self-consistent relativistic RPA 
approaches have been established by using the finite am- 
plitude method, where the i-FAM and m-FAM schemes 
are employed, respectively. 

For the FAM coding and calculations, the time-odd 
components of the functional, i.e., the nucleon currents 
and the space-component of the Coulomb field, must be 
kept explicitly. In the present covariant density func- 
tional, these time-odd components have the same cou- 
pling strengths as the corresponding time-even compo- 
nents due to the Lorentz symmetry. This makes the 
extension of FAM straightforward. Another key point 
for the FAM coding is the difference between the single- 
particle wave functions and their Hermitian conjugates. 
The formulas related to these key points are shown in 
Sec. In] in details. 

By taking the ISGMR in 208pb ^nd ^^^Sn as examples, 
the newly developed methods are verified by the conven- 
tional RPA calculations. It is also found that although 
the detailed shapes of the resonances depend on the box 



size R to some extents, the calculated centroid energies 
•nil /mo are precise up to 0.01 MeV. The experimental 
data in ^^^Pb is well reproduced. 

For the effects of the Dirac sea, it is confirmed that 
the ph configurations concerning the particle states in 
the Dirac sea must be included explicitly in the m-FAM 
scheme. On the other hand, such effects can be automati- 
cally taken into account in the coordinate-space represen- 
tation as in the i-FAM scheme, because the coordinate 
space, ^j. |r) (r| — ^ - \(j>j) {(j)j\, provides an equivalent 
complete set of basis for particle states. For the rear- 
rangement terms, instead of being calculated term by 
term in the conventional RPA, they can be implicitly cal- 
culated without extra computational costs in both i-FAM 
and m-FAM schemes. One simply needs to re-calculate 
the coupling strengths a and their derivatives da/dpb for 
each given set of (i/;'| and \i/j). 

In conclusion, the feasibility of the FAM for the co- 
variant density functionals has been demonstrated, and 
the advantages on treating the Dirac sea and rearrange- 
ment terms in the relativistic RPA have been presented. 
This opens a new door for developing the self-consistent 
relativistic RPA for deformed nuclei. 
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